Mode-coupling theory for reaction dynamics in liquids 
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A theory for chemical reaction dynamics in condensed phase systems based on the generalized 
Langevin formalism of Grote and Hynes is presented. A microscopic approach to calculate the 
dynamic friction is developed within the framework of a combination of kinetic and mode-coupling 
theories. The approach provides a powerful analytic tool to study chemical reactions in realistic 
condensed phase environments. The accuracy of the approach is tested for a model isomerization 
reaction in a Lennard- Jones fluid. Good agreement is obtained for the transmission coefficient 
at different solvent densities, in comparison with numerical simulations based on the reactive-flux 
approach. 

I. INTRODUCTION lates the diffusive Brownian motion of the reactive system 

near the barrier region. In this limit, the memory fric- 
tion is simply given by the Sjogren and Sjolander mode- 
coupling expression for the memory kernel of a GLE for 
self-motion in neat liquids 

In this paper we develop an alternative theory to 
treat non-Gaussian fluctuations within the GLE formal- 
equation./^ Grote and Hynes^ have extended this ap- ism. First, following the work of Oppenheim and his 



Chemical reactions in condensed phases are frequently 
described using two different methods^ The simple ap- 
proach consists of a Brownian particle moving in a one- 
dimensional bistable potential. In this case the dynam- 
ics of this system can be described by the Langevin 



proach to non-Markovian processes for a parabolic bar- 
rier, where the dynamics are given by the general- 
ized Langevin equation (GLE).^ These and other the- 
ories were tested numerically by Straub, Borkovec and 
BernCj^iSi and their results were crucial in the develop- 
ment of the PoUak, Grabert, Hanggi turnover theory^ 

The other approach, which became more common, is 
based on molecular dynamics simulations, and can ac- 
count for a general form of the Hamiltonian. Special 
methods have been developed to accelerate the barrier 
crossing and thereby numerically determine the rate con- 
stant. Chandler^ showed that in the time-correlation ap- 
proach to rate constants^Siifi the reactive flux rapidly de- 
cays to a plateau value, which can then be associated with 
the slow rate for crossing the barrier, '^-'^ For cases in which 
no preconceived notion of mechanism or transition state 
is known, a more general computational method, called 
the "transition path sampling," has recently been devel- 
oped and applied to complex reactive systems liSiiiiiii^ 

One of the major limitations of the simple approach 
based on the GLE formalism is related to how one can 
evaluate the frequency dependent memory friction (re- 
lated to the random force by the fluctuation dissipa- 
tion theorem). Most applications assume a Gaussian 
random force for the GLE. In this limit, the dynamics 
can be transformed into a Hamiltonian description where 
the system is linearly coupled to a harmonic bath^Siii 
Bagchi and his coworkers have used a mode-coupling the- 
Qj,yi8,i9,20 generalize this approach in order to include 
non-Gaussian fluctuationsiSiiS^ They assumed that the 
translational friction is the principal quantity that regu- 
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coworkers r^iSSiSS; we derive a GLE for the dynamics of a 
reactive system in a liquid host. Based on the formal ex- 
pression for the memory friction of this GLE, and using 
a combination of kinetic and mode-coupling theories, we 
obtain a simple expression for the memory friction. This 
memory friction is then used to obtain the reaction rate 
within the framework of the Grote- Hynes theory. Our ap- 
proach is different from that of Bagchi and coworkersSiiSS 
in two ways. First, we need not assume that the transla- 
tional friction is the principal quantity that regulates the 
diffusive Brownian motion of the reactive system. Sec- 
ond, it can be shown that this assumption limits the form 
of the couplings between the system and the liquid host. 
Thus, a more general form of the Hamiltonian can be 
treated within our formulation. 

Our paper is organized as follows: In Section^]we pro- 
vide a derivation of the GLE, an overview of the Grote- 
Hynes theory, and a derivation of our kinetic and mode- 
coupling theories. Tests of our approach for a model 
isomerization reaction in a LJ fluid is provided in Sec- 
tion lllll Finally, we conclude in Section Hvl 



II. THEORY 

The development of our approach to reaction dynam- 
ics is based on three steps. First, following the work of 
Oppenheim and his coworkers we derive a GLE 
for the dynamics of a subsystem. Second, we adopt the 
Grote-Hynes theory;^ which relates the reaction rate to 
the memory kernel of the GLE. Finally, we develop a the- 
ory based on a combination of kinetic and mode-coupling 
approachesi2ii2i2& to calculate the reaction rate within 
the Grote-Hynes formalism. 
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A. Generalized Langevin Equation 

Consider a general Hamiltonian for a reactive system 
in a liquid host of the from 



H = H,iR, P) + Hb{r, p) + ^{R, r), 



(1) 



where Hg {R, P) is the Hamiltonian of a reactive one di- 
mensional system with phase-space coordinates R and P, 
and reduced mass -ffb(r, p) is the Hamiltonian of N in- 
teracting liquid particles with phase-space coordinates r 
and p, and mass m, and (j){R^ r) is the coupling between 
the system and the solvent. 

Using the projection operator formalism of Zwanzig 
and MoripiSLSS and following the work of Mazur and 
Oppenheim(24 we now derive reduced equations of mo- 
tion for the system. First, we define a projection operator 
V that projects out the bath variables: 



PS = y drdpp{R,r,p)B = {B), 



(2) 



where 



p{R, r, p) = pb{r, p) exp{-/3[</)(i?, r) - w{R)]}. (3) 

In the above equation f3 — l/ksT is the inverse tem- 
perature, pb(r, p) = exp{— /?i/b(r, p)}/Z;, is the equihb- 
rium distribution function for the isolated bath, Zi, is the 
partition function of the isolated bath, and w{R) is the 
potential of mean force, given by: 

exp{-/3w(i?)} = J drdppb{r,p)exp{~l3(j){R,r)}. (4) 

Using the above projection operator it is straightfor- 
ward to show that the reduced equations of motion for 
the system are given by the GLE from:^i 



R 



P 

7 



(5) 



and 



P=-e 



\Lt 



dHs dw 



dR OR 



Kit) 



^ f dTe'^^P{KK{t-T)), 
M Jo 



(6) 



where, as noted above, p, is the reduced mass of the sys- 
tem, the (• • • ) is defined in Eq. jSJ, e*^* is the classical 
propagator with iL — {H, • • • }, and the random force 
K{t) is given by 



Kit) = 



d_ 
dR 



{cl>{R,r)-w{R)). 



(7) 



Note that the random force involves projected propaga- 
tion where the classical propagator e*^* is replaced with 
the projected propagator e*^^"^)^*. 



The memory kernel of the above GLE is related to the 
random force by the fluctuation dissipation theorem, and 
is given by 



C{t) = ^{KK{t)). 



(8) 



This correlation function is, in general, a function of time 
and a function of the system coordinate R (cf. Eq. (|7|)). 



B. Grote Hynes Theory 

To obtain the reaction rate one needs to solve the GLE 
given by Eq. However, as a result of the fact that 
the random force involves projected propagation of the 
couplings between the system and the bath, a complete 
solution of this GLE is an extremely difficult task. To 
circumvent this problem, Grote and Hynes have devel- 
oped an approximate theory to obtain the rate within 
the GLE formulationi^ 

The basic assumption made by Grote and Hynes is that 
in the barrier region, the dynamics can also be described 
by the above GLE, where the deterministic force (give 
[^ilt ~ f^]) replaced with an inverted parabolic 
approximation. After some lengthy algebra, Grote and 
Hynes obtained the following simple result for the trans- 
mission coefficient k: 



krsT 



(9) 



The transmission coefficient is given by the ratio of the 
reactive frequency (A^) to the barrier frequency (wf,). The 
former is related to the memory friction appearing in the 
GLE, and is given by 



(10) 



where C}{Xr) is the Laplace transform of C}{t) given by 



C\\r) = / dte~^^Wt), 



(11) 



and the symbol "f" indicates that the observable is eval- 
uated at the saddle point, namely, in the barrier region 
as required by the Grote-Hynes theory. To obtain the 
rate Eq. must be solved self-consistently. 

The Grote-Hynes expression for the rate reduces to 
Kramer's result if the memory friction decays rapidly 
on timescales faster than the system's motion.^ In the 
strong friction limit, Q{Xr)/p. » A,., the reactive fre- 
quency will be much smaller than the barrier frequency 
resulting in a small transmission coefficient. In the weak 
friction limit, Q{Xr)/p « A^, the reactive frequency 
simply equals the barrier frequency (A^ = w?,), and the 
rate is given by the transition state rate. This is a short- 
coming of the approximation made by Grote and Hynes 
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which fails to capture the turn over to Kramer's spatial 
diffusion regime. Hence, the Grote-Hynes rate expression 
is valid only for high frictions. 

The approach developed by Grotc and Hynes reduces 
the complexity of solving the GLE to that of estimating 
the memory friction of the GLE. For such problems, a 
very powerful approach has been developed and applied 
to many interesting dynamical problems in solutions. In 
the following subsection we describe a combination of 
kinetic and mode-coupling theories that we have devel- 
oped to model the memory friction of the GLE in order 
to study reaction dynamics in liquids. 



C. Kinetic and Mode-Coupling Theories 

One major difficulty in evaluating the memory ker- 
nel of the GLE given in Eq. ^ is that the cor- 
relation of the random force involves propagation of 
(f){R,r) — w{R) in the projected subspace spanned by 
Q = 1 — V (see Eq. ((TJ). To overcome this dif- 
ficulty we develop a theory to calculate C(i) that is 
based on a combination of kinetic and mode-coupling 
theoriesii^ This combination has been used to study den- 
sity and current fluctuations fiSiiSiSi solvation and relax- 
ation dynamics r2422i2£ and nonlinear spectroscopy2i*^SiS 
in classical liquids. A quantum mechanical generalization 
of this approach has recently been developed and applied 
successfully to study the dynamic response in quantum 
liquids |^4i2Si^Si2L^gi22ii2 

The Grote-Hynes approach requires that the memory 
friction be approximated in the barrier region. We there- 
fore replace the full memory friction of the GLE (C(0) 
with that approximated at the saddle point {(^t)), where 
(as before) the symbol "f" indicates that the position 
variable of the system is taken at the saddle point. Next, 
we replace (i) with an approximate form given by^ 

c\t)^cUt) + elicit) ^ (12) 

where Cb(^) aiid CMci^) are the "binary" and "mode- 
coupling" terms of the memory friction, respectively. The 
fast decaying binary term is determined from a short- 
time expansion (to second order in time) of the exact 
memory friction, and is given by 



and 



CL(t)-Ct(0) eM~{t/rn 



where the lifetime r is given by 



1 



2Ct(0)- 



(13) 



(14) 



In the above equations C^(0) and C(0) are the zero and 
second time moments of the memory friction, and are 
given by 



C^(o) 



f3 



d 



(15) 



Ct(0) = 



L /f 9 d 



\ dW dv 



(0(i?t,r)-«;(i?t)) 



(16) 



As noted above, is the reduced mass of the system, and 
m is the mass of a liquid particle. 

The slow decaying mode-coupling portion of the mem- 
ory kernel, Cmc(^)' niust be obtained from a mode- 
coupling approach. The basic idea behind this approach 
is that the random force projected correlation function 
decays at intermediate and long times predominantly into 
modes that are associated with quasi-conserved dynam- 
ical variables. It is reasonable to assume that the de- 
cay of the memory kernel at long times will be governed 
by those modes that have the longest relaxation time. 
In the present application the slow decay is basically 
attributed to couplings between wavevector-dependent 
density modes of the form 



where the self-density mode is given by 



Ck 



(17) 



(18) 



with Vj being the coordinate of liquid particle j , and the 
density mode is given by 



N 



(19) 



for N liquid particles. 

In practice, the simplest way to extract the dominant 
slow contribution of the decay of the memory friction is 
to introduce another projection operator, P2 given by 



5k.q 



kq 



NS{q) 



(&-k,-q, ■ 



(20) 



that projects any variables on the space spanned by 
the slow modes Ck and riq. Then, following the com- 
mon approximations in which the projected dynamics 
of the random force is replaced with the dynamics pro- 
jected onto these slow variables, and replacing four-point 
density correlations with a product of two-point density 
correlationSfi2ii2iS& we find that Cmc(^) given by 



JL 
lip 



(27r)3 
F,{q,t)F{q,t) - 



|yt(q)|V5(<z)2] 
' F,biq,t)Fh{q,t)] 



(21) 



In the above equation p is the liquid number density, 
S{q) is the structure factor of the neat fluid, Fs{q,t) and 
F{q, t) are the self-intermediate and intermediate scat- 
tering functions of the neat fluid, respectively, and the 
vertex is given by the static average of the product of 
the coupling force and the slow modes: 



^^(q) = (^[0(i?^r)-^i?t)]5^_^ 



(22) 
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As before, the symbol "f" indicates that the couphng 
force is evaluated at the saddle point. 

The binary self-intermediate and binary intermediate 
scattering functions are given by 



The solvent Hamiltonian is given by the LJ form: 



and 



Fsbiq,t) = exp 



Fh(g,t) = 5(g)exp 



2(3mS{q) j 



(24) 



The subtraction of the product of these terms in Eq. 1)2 l|l 
is done to prevent over-counting the total memory kernel 
at short times, namely, to ensures that the even time 
moments of the total memory kernel are exact to forth 
order in time. 

To obtain the memory friction one requires as input the 
first two moments C^(0) and C^(0): the static structure 
factor S{q) and the vertex T^^(q). These static averages 
can be obtained from simulations, or from the proper 
integral equation formulation. In addition, we need the 
self-intermediate and intermediate scattering functions. 
These time-dependent correlation functions of the neat 
fluid can be obtained from simulations, or alternatively 
from a similar mode-coupling approach, where a GLE for 
the density fluctuations is solved using a combination of 
kinetic and mode-coupling theories>i2ii2iS£ 



III. RESULTS 

To assess the accuracy of the proposed theory we have 
studied the rate of a model for an isomerization reaction 
of a diatomic molecule in a Lennard- Jones (LJ) fluid. 
The Hamiltonian of the entire system and bath can be 
described by Eq. Q . The isomerizing diatomic molecule 
is made of two atoms with equal mass m*, interacting 
via a symmetric double-well potential. Without loss of 
generality, we place the atoms along the 2;-axis, where the 
position vectors are Ri = {0,0,0} and R2 = {0,0, i?} 
for atom 1 and 2, respectively. The reaction coordinate 
is taken as the distance R = |Ri — R2I separating the 
atoms. We allow the atoms to move only along the z axis 
(which is the reaction coordinate). The reactive system 
Hamiltonian can be described by: 



Hs{R, P) 



p2 

2^ 



+ V^{A[{R-R^)/af-l}\ (25) 



where = m*/2 is the reduced mass of the system, Vq 
is energy barrier separating reactants from products, R^ 
is the location of the transition state, and a is the dis- 
tance between the two minima corresponding to stable 
reactants and products. For the results shown below we 
take i?t = 21/6(7, Vq = bksT, T = 2.5e, a = 2a/3, and 
a reduced mass fJ- = 1/2 for the isomerizing diatomic 
molecule. 



N 



^-^ 2m 



N 

E 

i>j=i 



4e 



(26) 

(23) where Vj is the position vector of liquid particle j with 
momentum pj and mass m, and — \vj — r^j. 

For simplicity we take the solute-solvent interaction 
to be a two-site LJ potential in which each site of the 
molecule interacts with the solvent atoms through pre- 
cisely the same LJ potential: 



(R,r) 



N 

E 



4e 



E 



4e 



- 



a 

"23 



(27) 



Here Vij — \Yj 



Ri|. In the application of the mode- 
coupling theory we need to evaluate the derivative of 
(/)(R, r) with respect to the reaction coordinate R. This 
is done using the chain rule, and the results is given by 



90(R,r) 
dR 



1 p0(R,r) 

2 1 dZ2 



90(R,r) 
dZi 



(28) 



where Zi and Z2 are the z components of the position 
vector of the diatomic molecule Ri and R2, respectively, 
and the factor i comes from the Jacobian. The deriva- 
tives of 0(R, r) along the other directions vanish since 
the system is frozen along the z-axis. 

To test our approach, we have calculated the static and 
dynamic input required to obtain the time-dependent 
friction, using the Monte Carlo and molecular dynamics 
simulation techniques. While other approaches, such as 
the integral equation theory and a proper mode-coupling 
treatment of density fluctuations, can be used to obtain 
the static and dynamic input, we feel that a fare test of 
our theory should rely on numerically exact input. Thus, 
despite the success of the theoretical approach in predict- 
ing structural**^ and dynamical^^ properties in LJ sys- 
tems, we limit our study to the more accurate simulation 
approach. 

The static input was obtained using the molecular 
dynamics (MD) method for — 500 particles (in- 
cluding the diatomic) and for a set of densities p = 
0.5,0.6,0.7,0.8,0.9,1.0 in reduced LJ units# Approxi- 
mately 10*^ MD steps were made for each density. Every 
50 steps we collected data for C^(0), S{q) and for 

the vertex V'^ (q) . The self-intermediate and intermediate 
scattering functions were calculated using the molecular 
dynamics method for identical conditions. The results 
for each density were averaged over 10 different runs of 
total run time of < = 15 reduced LJ units. 

In Fig. n] we plot the static structure factor and 
the vertex squared for all densities studied in this 
work. The features observed for the static structure 
factor are well understood, and have been discussed 
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FIG. 1: The vertex square (left panel) and the static struc- 
ture factor (right panel) for different liquid densities p = 
0.5, 0.6, 0.7, 0.8, 0.9, 1.0 at T = 2.5. All results are in reduced 
LJ units. For clarity, the results for different densities are 
shifted vertically, from low to high densities. 



elsewhere 1^ The magnitude of the vertex given by 
l/t(q) = (g|^[(/,(Rt,r) -u;(Rt)]6q,-q) determines the 
contribution of the different density modes to the decay 
of the memory friction at intermediate and long times. 
We find that at high liquid densities the major contribu- 
tion comes from modes with a characteristic wavelength 
oi q K, 27r/(T. This wavelength corresponds to the average 
intcrparticlc separation, and its value slightly increases 
with decreasing density (note the shift in the position 
of the maximum of S{q) with density). Thus, density 
modes associated with liquid motion on length scales of 
the interparticle separation (namely, on length scale of a) 
contribute the most to the friction at intermediate and 
long times. 

As expected, the contribution of lower wavelengths to 
the decay of the memory friction becomes more signif- 
icant at lower liquid densities. In fact, at the lowest 
density studied, the contribution of density modes be- 
low q K, lix I a to the decay of the memory friction at 
intermediate and long times is more significant than the 
contribution of modes near the first peak in Sici). We 
attribute this effect to the change in the mechanism of 
self-diffusion of liquid particles. At high liquid densi- 
ties, the self-diffusion is dominated by opening of the 
cage surrounding a system, with a typical length scale 
of cr related to the size of the cage. This is not true at 
low liquid densities, where liquid particles can hop over 
much larger distances. In addition, we find that motion 
within the cage become significant at higher liquid den- 
sities, as refiected in non vanishing values of the vertex 
above qK,2-n ja. 

Using the static input obtained from the Monte Carlo 
simulations along with the values of the self-intermediate 
and intermediate scattering functions obtained from 
molecular dynamics simulations, we have calculated the 




time 



FIG. 2: A plot of the normalized memory friction (t) 
given by Eqs. and H21^ for different densities p = 

0.5, 0.6, 0.7, 0.8, 0.9, 1.0 at T = 2.5, in reduced LJ units. 



total memory friction given by Eqs. 1)13(1 and \2,\\ . The 
results are shown in Fig. (j^J for all densities studied in 
this work. As can be seen, the decay of the memory 
friction is characterized by two time scales. A fast decay 
dominated by C}^ {€) followed by a slower decay dominated 

by Cmc(^)- "^^"3 contribution of the mode-coupling por- 
tion to the total memory friction is significant only at 
high liquid densities. At the lower densities studied the 
mode-coupling portion of the memory friction is negligi- 
ble. 

For the specific model studied here where the solute- 
solvent interactions equal to the solvent-solvent interac- 
tions, the memory friction can also be obtained by invert- 
ing the GLE for the velocity autocorrelation function of 
a neat fluid/dS. Straub, Borkovec, and Berne^^ have calcu- 
lated the velocity autocorrelation function for the LJ fluid 
using the molecular dynamics technique, and obtained 
the memory friction by inverting the proper GLE. Our 
results obtained from the mode-coupling theory agree 
well with their simulated results (not shown). Specifi- 
cally, we find that the decay of the memory friction at 
short times is nearly independent on the liquid density, in 
agreement with the molecular dynamics resultst^ More- 
over, at high liquid densities, our theory provides semi- 
quantitative agreement with the molecular dynamics re- 
sults at all times. For low liquid densities, we observe 
small deviations from the simulation results at interme- 
diate times. This shortcoming of the mode-coupling ap- 
proach is expected since the simulated memory friction 
becomes slightly negative at intermediate times, and the 
mode-coupling approximation is known to fail under such 
circumstances. However, since the contribution of the 
mode-coupling portion to the memory friction is rela- 
tively small at these low densities, this has a vanishing 
effect on the value of the transmission coefficient. 
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t(0) 



FIG. 3: The transmission coefficient as a function of the in- 
tegrated friction in the limit of spatial diffusion. The stars 
connected by solid line are the results of the present mode- 
coupling theory. Solid line is the result of the weak collision 
theory using an exponential frictiouf^ and the dashed line is 
the Kramer's result. 

Using the memory friction obtained from the kinetic 
and mode-coupling theory we have calculated the trans- 
mission coefficient by solving Eq. Q self-consistently. 
The results for Kj 0jS Oil function of the integrated fric- 
tion Ct(0) (cf. Eq. ((TTJ) are shown in Fig. O We com- 
pare our results to the results obtained by the weak colli- 
sion theory (based on a simple connection formula) using 
an exponential friction,''^ and to Kramer's theory.^ For 
the present model the weak collision theory provides a 
quantitative agreement for the transmission coefficient 
(within the noise level of the simulations) in compari- 
son with results obtained using the absorbing boundary 
approximation^2i^ to the reactive-flux formalism iSiiS As 
can be seen, our theory provides quantitative results for 
the transmission coefficient over the entire range of fric- 
tions studied. 

The different values of the integrated friction CU^) at 
which we have calculated k, where obtained by scanning 
the liquid density from low p = 0.5 to high p = 1.0 
values. For each density we have calculated the static 
and dynamic input required by our mode-coupling the- 
ory, and obtained the transmission coefficient by solving 
Eq. self-consistently. Alternatively, one can control 
the integrated friction (0) by changing the value of the 
reduced mass /i of the isomerizing diatomic molecule (cf. 



Eq. O. We find that our mode-coupling theory provides 
similar quantitative results (not shown) for the transmis- 
sion coefficient when the reduced mass of the isomerizing 
diatomic molecule is varied. 

IV. CONCLUSIONS 

We have presented a theoretical approach for the cal- 
culation of reaction rates in condensed phases based on 
the Grote-Hynes formalism. A combination of kinetic 
and mode-coupling theories were developed to obtain 
the memory friction required by the Grote-Hynes for- 
malism. The approach was applied to study a model 
isomerization reaction of a diatomic molecule in a LJ 
fluid. Good agreement for the transmission coefficient 
was obtained in comparison with the simulation results 
of Straub, Borkovec, and Berne^ based on the reactive 
flux formalism. 

Unlike simulation techniques, our approach is a theory 
and thus provides additional insight into the reaction dy- 
namics in liquids. For example, we showed that the con- 
tribution of the mode-coupling portion to the decay of the 
memory kernel at intermediate times is significant only 
at high liquid densities. Thus, an accurate description 
of the reaction dynamics at low liquid densities can be 
obtained from a kinetic theory alone. Furthermore, the 
mechanism for the decay of the memory friction (which 
is reflected in the reaction rate) is quite different at low 
versus high liquid densities. At high liquid densities the 
decay of {t) is dominated by liquid modes with a length 
scale comparable to the separation between the fluid par- 
ticles, while at lower liquid densities motion on larger 
length scales also contributes the decay of C^(^)- This is 
significant for the development of coarse grained models 
for reaction dynamics. 

We believe that our approach will be useful in other sit- 
uations in which simulation techniques are still limited. 
For example, for reaction dynamics in supercooled liq- 
uids that are characterized by slow density fluctuations. 
Or for liquid hosts that are characterized by quantum 
mechanical susceptibilities. Work along these directions 
is currently underway. 
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